/*******************************************************************************

                                BACK

  Title:        Back transformation of calculated diffraction patterns.    
  Author:       Hanna Szoke
  Date:         11/30/92
  Function:     This program estimates electron densities from a set of
		calculated diffraction patterns with phases. It obtains
		a "solution map" - i.e., the amplitudes of a set of
		Gaussian densities of given width, eta*(grid_spacing)^2,
		centered on a simple grid of given grid spacing. 

		It is assumed that the input fc file has been correctly
		apodized prior to running Back.

		One purpose of this calculation is to provide Solve 
		with a map of initial values for the conjugate gradient 
		optimization process which is everywhere non-negative
		and consistent with the reciprocal-space representation.

		Output includes a log of the process, together with a
		binary file of electrons/voxel.

*******************************************************************************/

#include "speden.h"	/* includes util.h, cellparams.h, symmetry.h, dims.h */
#include "speden_Np.h"
#include "cost.h"
#include "ccg.h"		/* needed for dfdxpc */
#include "sp.h"

COMPLEX	*Fc0 ;		/* Calculated structure factor for known atoms */

void    convolve(int, real *, real *, COMPLEX *) ;

static	char	b_title[MAXSTRING] ;	/* whatever the user wants ... */

static	void	back_setup() ;
static	void	b_alloc_hklArrays() ;
static	void	report_bmemory() ;
static	void	b_readInput() ;
static	void	b_echo() ;


void	back_main(int argc, char *argv[])
{
     void	minimal_Lb_funct(real *) ;	
     void	do_inplace_cs_symmetrization(real *, real *) ;
     void	Lb_get1Rfac(COMPLEX *) ;	
     void	shellHeader_R() ;				
     void      	b_HoloMethod() ;		/* in ccgutil.c */
     void	analyzEpixFile(real *) ;	/* in qshare.c */
     void	writeEpixFile(char *, real *) ;	/* in qshare.c */
     void	write_HDFnump(char *, real *) ;
     void	write_HDFmap(char *, real *, real, real) ;
     void	transform() ;
     int        get_min_max(real *, int, real *, real *) ;
     void	printXmap() ;

     double	sum = 0 ;		/* NOTE: 16/2/2005	*/
     real	residue = 0 ;			
     int        n ;
     char	problemName[MAXSTRING] = "" ;

  /******************************************************

  Check for right ballpark; 2/25/02: sf_filename may be 
  entered EITHER on execute line OR from input file!
  ******************************************************/

     strcpy(sf_filename, "none") ;

     if (argc > optind+3)  	/* too many arguments */
          ballpark(caller) ;

     if (argc > optind+2)
	  sprintf(sf_filename, argv[optind+2]) ;

  /*********************************************
  Pull off prefix defining whereabouts of input 
  file and remove .inp extension.
  Say hello, prepare cost file.                             
  *********************************************/

     strip_path(input_filename, pwd_name) ;
     strncpy(problemName, pwd_name, strlen(pwd_name) - 4) ;

     hello(caller) ;

     if (verbose)
          costfile_option = TRUE ;
     else
          costfile_option = FALSE ;

     if (costfile_option) {
	  sprintf(cost_filename, "%s_back.cost", problemName) ;

          if ((fp_cost = fopen(cost_filename,"w")) == NULL) 
          {
               sprintf(message, "Cannot open %s", cost_filename);
               Fatal(message) ;
          }
     }

  /*********************************************
  Fetch input conditions, do all initializations.
  Report initial standard deviation 
  *********************************************/

     back_setup() ;
     minimal_Lb_funct(&residue) ;

     sprintf(message, 
	 "\nInitial standard deviation = %g", sqrt(residue / (Nfc-1))) ;
     
     printTwice(message) ;

  /*********************************************
  Do holographic reconstruction.                           
  *********************************************/

     sprintf(message, "\nApplying phased holographic reconstruction \n") ;
     printTwice(message) ;

     b_HoloMethod();

  /*********************************************
  Symmetrize resulting Np arrays.                         
  Report final standard deviation.
  *********************************************/

     if (Ncs > 1)  
          do_inplace_cs_symmetrization(nump, nump) ; 

     minimal_Lb_funct(&residue) ;

  /*********************************************
  Summarize results, write electron/pixel arrays.
  *********************************************/

     for (n = 0; n < Nptotal; n++) 
          sum += *(nump + n) ;
          
     sprintf(message, "Sum of recovered electrons is %g\n", sum) ;
     printTwice(message) ;

     analyzEpixFile(nump) ;
     convolve(Nfc, nump, expfac, Oar);
     shellHeader_R() ;
     Lb_get1Rfac(Oar) ; 
     e_free(Nfc*sizeof(COMPLEX), Oar) ;

     sprintf(out_filename, "%s_back", problemName) ;
     writeEpixFile(out_filename, nump) ; 

  /*********************************************
  Transform (in-place) .bin info to map.
  Note: in addition to the regular map (in xplor 
  or hdf format), we also write an ordinary .bin
  file for other post-processors.
  *********************************************/

     transform(Nptotal, nump, 0) ;
     get_min_max(nump, Nptotal, &ar_min, &ar_max) ;

     if (HDFOUT) {
          write_HDFmap(out_filename, nump, ar_min, ar_max) ;
     }
     else {
          sprintf(message," %7d %7d %7d %7d %7d %7d %7d %7d %7d",
             Nx, 0, Nx-1, Ny, 0, Ny-1, Nz, 0, Nz-1) ;
          printXmap(out_filename, message, Nx, Ny, Nz, nump) ;
     }
     strcat(out_filename,"_map") ;
     writeEpixFile(out_filename, nump) ;

     return ;
}

void	back_setup()
{
     void	setup_csind() ;			/* in crysutil.c */
     void	hklTabHead() ;			/* in hklutil.c */
     void	fsetDsqLimit(int, char) ;	/* in hklutil.c */
     int	readFcList(char *) ;		/* in hklread.c */
     int	read_HDFfcalc(char *) ;		/* in hdfutil.c	*/
     real	*setupexpfac(real) ;		/* in hklutil.c */
     void	echo_ignored() ;		/* in init.c	*/
     void	initMinMaxType(int) ;		/* in ccgutil.c */
     void	alloc_ccg_mem(int) ;		/* in ccg.c */
     void	initfft(int, int) ; 		/* in fftshare.c */
     void	readEpixFile(char *, real *) ; 

     int        n ;

  /*********************************************
  Read input parameter file.
  *********************************************/

     b_readInput() ;
     b_echo() ;
     echo_ignored() ;
     report_bmemory() ;

     printTwice("Setting up initial arrays ...") ;

     b_alloc_hklArrays() ;

  /*********************************************
  Read a set of diffraction patterns, Fc0.
  Prepare exp factors for FFTs using delfac = 
  PISQ*eta*dr*dr.
  *********************************************/

     if (HDFIN)
          Nfc = read_HDFfcalc(sf_filename) ;
     else
          Nfc = readFcList(sf_filename) ;

     initfft(1, Nfc) ;

     expfac = setupexpfac(delfac) ;

     fsetDsqLimit(Nfc, FALSE) ;
     hklTabHead() ;

  /*********************************************
  Identify scope of problem.                               
  *********************************************/

     sprintf(message,
           "\n\t\tThis problem has %d equations, %d unknowns\n", Nfc, Nptotal) ;
     printTwice(message) ;

  /*********************************************
  Determine whether or not there is a starting 
  model and, if there is, read it into the knownp.
  Prepare arrays for the electron densities and
  for the conjugate gradient solver.  
  *********************************************/

     nump = (real *) e_malloc(Nptotal*sizeof(real), "back_setup") ;

     if (strcmp(md_filename,"empty") != 0) {
         knownp = (real *) e_malloc(Nptotal*sizeof(real), "back_setup") ;
         readEpixFile(md_filename, knownp) ;
         for (n = 0; n < Nptotal; n++)
	     *(nump+n) = *(knownp+n) ;
     }
     else 
         for (n = 0; n < Nptotal; n++)
	     *(nump+n) = 0 ;
     
     initMinMaxType(Nptotal) ;
     alloc_ccg_mem(Nptotal) ;

     if (Ncs > 1)  
          setup_csind() ;
     
     return ;
}

void	b_alloc_hklArrays()
{
     void       alloc_fc(int) ;

  /*******************************************************
  Set up arrays in (h,k,l) reciprocal lattice triplets.
  Prepare for FFTs.  Note that the initial Nfc (length of 
  file) is corrected DOWN at read time because some input 
  may not actually be used.  
  *******************************************************/

     alloc_fc(Nfc) ;

     strcpy(message, "b_alloc_hklArrays") ;
     Oar   = (COMPLEX *) e_malloc(Nfc*sizeof(COMPLEX), message) ;
}

void	report_bmemory()
{
     int	HDFdatalength(char *, char *) ;
     real	mem_Npspace, mem_recspace, mem_tot ;

     printTwice("Checking out input fcalc file\n") ;

     if (HDFIN)
          Nfc = HDFdatalength(sf_filename, "fcalc/Nfc") ;
     else
          Nfc = filelength(sf_filename) ;

	/*******************************************************
	In Np space, there are real arrays (nump, minp and maxp),
	1 int array (typep),possibly a knownp, and 
	4 CCG real arrays (xn, gn, go and d), each with Nptotal 
	elements.  There is also a complex FFT array in Nptotal.
	*******************************************************/
     
     mem_Npspace = (8*sizeof(real) + sizeof(int) 
                    + sizeof(COMPLEX)) * Nptotal ;

	/*******************************************************
	In reciprocal space, there are 3 complex arrays in 
        Nfc (Fc0, Oar, and f for NDFT) and 7 real arrays (h,k,l 
        lists as read in and for NDFT and expfac).  
	*******************************************************/

     mem_recspace = 3 * Nfc * sizeof(COMPLEX) +
		    7 * Nfc * sizeof(real) ; 

     mem_tot = ((mem_Npspace + mem_recspace) * 1e-9) + 0.5 ;
     if (mem_tot > 1.) 
           sprintf(message, 
     "Approximate memory requirements in Gbytes: %6.2g", mem_tot) ;

     else {
          mem_tot = ((mem_Npspace + mem_recspace) *1e-6) + 0.5 ;
          sprintf(message, 
     "Approximate memory requirements in Mbytes: %6.2g", mem_tot) ;
     }
     printTwice(message) ;
}


#define	DFDX_BACK	1.e-2    /* changed from 1.e-3, June 30, 2005 */

void	b_readInput()    
{
     int	k ;

     /************************************************************
     Read basics; set defaults for global variables.
     ************************************************************/

     readBasicInput() ;

     dfdxpc = DFDX_BACK ;
     discrp_frac = 1. ;

     strcpy(b_title, "") ;
     strcpy(md_filename, "empty") ;

     /************************************************************
     Now read input that is specific to Back
     *************************************************************/

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "TITLE") == 0) {
               strcpy(b_title, nextline+6) ;
	       goodline[k] = TRUE ;
          }
          else if ((strcmp(id, "FC_FILENAME") == 0) ||
              (strcmp(id, "FC_FN") == 0)) {

               if (strcmp(sf_filename, "none") == 0) {
                    sscanf(nextline, "%*s %s", sf_filename) ;
	            goodline[k] = TRUE ;
	       }
          }
          else if ((strcmp(id, "MD_FILENAME") == 0) ||
                   (strcmp(id, "MD_FN")       == 0)) {
               sscanf(nextline, "%*s %s", md_filename) ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "DFDX_CRIT") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               dfdxpc = t1 ;
	       goodline[k] = TRUE ;
          }
     }

     /************************************************************
     Check for fatally flawed input.
     ************************************************************/

     if (strcmp(sf_filename, "none") == 0)
          Fatal("Missing fc file name!") ;

}

				/*********************************************
          			Echo information in .inp file to log.     
				*********************************************/
void b_echo()      
{
     if ((int)strlen(b_title) > 0) 
          fprintf(fp_log, "\n%s\n\n", b_title) ;
     else
          fprintf(fp_log, "\n") ;
     
     if ((int)strlen(sf_filename) < 40)
          fprintf(fp_log, 
          "Calculated structure factors will be read from %s\n", sf_filename) ;
     else {
          fprintf(fp_log, "Calculated structure factors will be read from\n") ;
          fprintf(fp_log, " %s\n", sf_filename) ;
     }
     if (strcmp(md_filename,"empty") != 0) {
          if ((int)strlen(md_filename) < 40)
               fprintf(fp_log, 
     "A starting physical space model will be read from %s\n", md_filename) ;
          else 
               fprintf(fp_log, 
     "A starting physical space model will be read from\n %s\n", md_filename) ;
     }
     else
          fprintf(fp_log, "There is no starting physical model\n") ;

     fprintf(fp_log, "\n") ;
     fflush(fp_log) ;
}
